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.1.  Introduction  and  Objectives; 
1 . 1  introduction s 


j1\o  nonlinear  analysis  of  roinfom'd  a  morel  o  MtnwluivM  by  the  T  ini.Uc^ 
element  method  cannot  be  succor  fully  performed  if  the  principal  sources  of  ma¬ 
terial  nonlinear  behavior  are  not  included  in  tlie  formulation.  The  material 
characteristics  that  have  to  be  considered  are  the  nonlinear  stress  strain  re¬ 
lation  Cor  the  concrete,  the  stress  strain  relation  for  Uk:  reinforcement, 
concrete  anisotropy  due  to  complex  stress  states  and  cracking ,  postcracking 
shear  transfer  mechanisms  at  open  cracks,  and  the  concrete  reinforcement  bond 
slip  relations. 

This  report  presents  a  conputer  subroutine  for  the  nonlinear  analysis  of 
reinforced  concrete  elements  that  includes  the  nonlinear  stress  strain  rela¬ 
tions  for  the  concrete  and  the  reinforcement,  concrete  anisotropy  due  to 
cracking  and  multiaxial  stress  states,  and  the  postcracking  shear  transfer 
mechanisms  present  at  a  slightly  open  crack.  The  nonlinear  behavior  of  the 
concrete  is  represented  by  the  endochronic  model  (  5  )  while  the  stress  strain 
relation  for  the  reinforcement  represents  the  elastic,  plastic,  and  strain 
hardening  stages  under  monotonic  or  repeated  loads.  The  postcracking  shear 
transfer  mechanisms  included  in  tlie  subroutine  are  the  interface  shear  trans¬ 
fer  mechanism  on  the  rough  surfaces  of  a  cracked  plane  and  the  dowel  action  of 
the  reinforcement  crossing  the  crack.  The  subroutine  developed  uses  the 
distributed  crack  approach  to  combine  the  stiffness  matrix  for  the  uncracked 
concrete  with  the  stiffness  relation  for  the  cracks  that  predicts  the  incremen¬ 
tal  stresses  induced  in  a  reinforced  concrete  finite  element  by  a  set  of 
prescribed  incremental  strains. 

Ti«  different  sections  of  this  report  describe  first  the  constitutive 
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a  detailed  discussion  of  Ml  the  axle  subroutines.  Finally,  several  experimen¬ 
tal  tests  are  compared  with  the  results  predicted  1  jy  the  comjxjtor  prtxjram  to 
determine  its  validity. 

1.2  Objectives: 

The  principal  objective  of  this  project  was  to  develop  a  material 
subroutine  that  included  the  princi^l  sources  of  nonlinear  behavior  in  rein¬ 
forced  concrete.  Hie  specific  objectives  arc: 

A.  Development  of  a  material  subroutine  that  calculates  the  incremental 
stress  vector  caused  by  a  given  vector  of  strains  in  a  plane  stress,  strain 
or  axisymetric  finite  element.  The  subroutine  should  consider  the  nonlinear 
behavior  caused  by  the  following  sources: 

1.  Stress  strain  relation  for  concrete  based  on  the  endochronic 
model  presented  by  Bazant  (2,3  ) . 

2.  Concrete  anisotropy  caused  by  cracking  and  multiaxial  stress 

states * 

3.  Fostcracking  shear  transfer  mechanisms.  Both  the  interface 
shear  transfer  and  the  dowel  action  stiffness  representation  are  ihcluded  in 
the  subroutine. 

4.  Stress-strain  relation  for  reinforcement  that  models  the  elastic, 
plastic,  and  strain  hardening  stages  for  monotonic  and  repeated  loads. 

B.  Comparison  of  the  material  subroutine  code  predictions  with  available 
experimental  data  to  establish  the  valitidity  of  the  proposed  formulation. 
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2 -  Constitutive  Relations  for  Nonlinear  Analysis  of  Reinforced  Concrete 
2 . 1  Introduction : 


The  following  sections  present  the  theoretical  background  required  to  model 
the  nonlinear  behavior  of  reinforced  concrete  in  the  computer  program.  First, 
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miW,  are  discunwil,  tolh'wixl  l»y  IN'  i-onstHutivo  mini  ions  t'n^ 'loyixl  U>  mnk'1 
tins  Ixilviviur  of  the  rt?infot omik'mI.  sub jested  to  mmotnnie  or  roproaU'd  loads. 

|'oi‘  oUmonts  that  turn?  not  maekid,  t  hr  coiuMvto  and  alivl  st  if fiK'Ss  tvlat ions 
In  tonus  of  tlte  glnlwil  nundinal  on  ran  l*'  iKVkil  to  obtain  an  ineroni'nl  al  si  ress- 
strain  relation  for  the  reinforced  concrete  element.  For  elements  that  have 
cracked,  however,  the  incremental  stress-strain  relation  has  to  consider  the 
constitutive  relation  for  ttie  cracks,  presented  on  the  last  section  of  this 
chapter. 

2 . 2  Constitutive  Relations  for  Concrete: 

Several  theories  have  been  dev loped  to  predict  the  response  of  plain  con¬ 
crete  to  multiaxial  stress  states  among  which  are  ttve  linear  and  nonlinear 
elasticity  theories,  the* work  hardening  plasticity  theories,  the  plastic  frac¬ 
turing  theory,  and  the  endochronic  theory  (4,8).  Of  tnese  theories,  the 
endochronic  theory  has  received  particular  attention  as  it  provides  a  contin- 
ous  model  for  the  nonlinear  representation  of  concrete  without  the  explicit 
formulation  of  a  yiel*  condition  and  hardening  rules.  The  endochronic  model 
developed  by  Bazant  and  co-workers  (2,5  )  have  been  used  succesfully  to 
predict  the  nonlinear  stress*  strain  curve  for  concrete  subjected  to  monotonic 
or  repeated  loading. 

The  endochronic  dieory  for  concrete  initially  proposed  by  Bazant  (  2  ) 
introduced  a  non  decreasing  scalar  variable,  denominated  intrinsic  time,  to 
represent  the  accumulation  of  inelastic  strains  as  a  function  of  the  strain 
increments  applied  to  the  element.  The  intrinsitive  increments  were  assumed 
to  be  sensitive  to  tnc  hydrostatic  pressure.  .  The  theory  also  modelled  the 
strain  hardening  and  strain  softening  regions  of  the  stress  strain  curve  for 
concrete,  the  inelastic  dilatancy  due  to  shear  straining  measured  by  another 
non-decreasing  scalar  variable,  and  the  dopendanee  of  the  incremental 
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oWmtic  ntxlull  on  Mk»  dUnl-mcy  iw'uunrc. 

'Hkj  initial  endochronic  mrxlel  pvtT}xvnod  Ijy  llixant  was  Inter  refined  (5  ) , 
to  include  the  additional  inelastic  strains  caused  by  hydrostatic  compression, 
while  changeR  in  tlie  strain  softening  ramie  of  the  stress  strain  curv^, 
de| Tendance  of  material  jvmimntors  (Xi  strength,  and  an  inpmvixl  description  of 
the  strain  softening  behavior  under  monotic  or  repeated  loading.  The  refined 
endochronic  model  was  used  to  represent  the  nonlinear  behavior  of  concrete  in 
the  material  subroutine . 


2.2.1  Endochronic  Model  for  Trialxial  Behavior  of  Concrete 


Hie  stress  strain  relations  for  the  endochronic  model  are  qiven  in  terms 
of  the  deviatoric  and  volumetric  relations,  as  follows: 


1  I 


Aeij  =  ASi^  +  Ae^ 


At: 


»  Ao  +  At 

3K 


Ac,  Ao 


Where:  Aei3,  AS^j  =  deviatoric  components  of  strain  and  stress  ■  tensor , 

respectively. 

*  volumetric  component  of  strain  and  stress  tensor, 
respectively. 

=  inelastic  deviator  strain  increment. 

=  inelastic  volumetric  strain  increment 
=  Bulk  and  shear  modulus. 

=  Cartesian  coordinates  indexes. 
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Tbe  volumetric  cortpononts  of  the  strain  and  stress  vectors  are  computed 


A  r. 


F11  +  F22  +  f33 
Ao  =  ou  +  o22  +  o33 

while  the  deviatoric  ccuponents  are  obtairied  from: 


2a 


2b 


The  inelastic  volumetric  strain  increment  is  a  function  of  the  volumetric 
stress,  o,  the  inelastic  dilatancy  X,  the  shear  compaction  X* ,  and  of  the  com¬ 
paction  intrinsic  time  parameter  Z* .  The  compaction  intrinsic  time  parameter 
z’have  been  introduced  to  account  for  the  volumetric  inelastic  strains  caused 
by  hydrostatic  stress  states,  while  the  shear  compaction  parameter  X*  accounts 
for  the  increased  volumetric  strains  observed  in  triaxial  stress  .tests  when 
compared  to  hydrostatic  stress  tests.  The  inelastic  voluretric  strain  is  given 
by: 

Ac”  =  AX  +  A  Z’  +  AX’  '  5 

Thus,  the  inelastic  deviator ic  and  volumetric  strains  are  a  function  of  the 
distortion  intrinsic  time  Z ,  the  compaction  intrinsic  time  Z  * ,  the  inelastic 
dilatancy  X,  shear  compaction  X1,  the  bulk  and  shear  modulus,  and  the  volume¬ 
tric  and  deviatoric  stress  components  present  in  the  element.  The  endochronic 
parameters  and  the  bulk  and  shear  modulus  are  computed  from  the  set  of  equa¬ 
tions  summarized  in  Appendix  Al.  It  should  be  noted  that  the  functions  used  to 
calculate  the  intrinsic  time  parameters,  inelastic  dilatancy  and  shear  compact¬ 
ion  are  a  function  of  the  current  stress  and  strain  invariants,  and  of  the 
principal  stresses  in  the  element.  Therefore,  for  a  prescribed  strain 
increment#  the  associatted  stress  increment  has  to  be  computed  in  an  iterative 


fashion  until  llw  « *iu  |»  m  •!  u « »t » i  • '  |  1 1 .  hi*  >  I  < -i  ■:  ,it  .my  i|iven  iloi.iliun  *  •»  »nvi  •* « §»  •  hi 
'thtWO  of  tile'  previous  floral  ion.  In  LIh'  computer  code,  Uio  iteration  is  termi¬ 
nated  when  the;  difference  lx? tween  the  previous  and  current  values  of  the  inelas¬ 
tic  dilatancy  and  shear  compaction  parameters  in  within  0.01. 

■  Tt  nliould  lx?  noUsI  that  the  hceuimilnlod  values  of  several  eudoehronic 
parameters,  as  well  as  for  the  stress  and  strain  vectors,  are  required  for  the 
equations  given  in  Appendix  Al.  The  current  value  of  any  parameter  at , the  end 
of  a  strain  increment  is  calculated  from  the  increment  of  said  parameter  computed 
when  the  iteration  is  finished,  and  the  value  obtained  in  the  previous  strain 
increment. 

The  stress  strain  relations  in  terms  of  deviatoric  and  volumetric  compo¬ 
nents  given  by  Equation  1  can  be  combined  to  obtain  an  incremental  stress  strain 
relation  in  terms  of  the  element  coordinates,  using  the  relations  given  by 
Equation  3.  Rearranging  Equations  1,  we  have 

I  I 

AS^j  =  2G  Aejj  —  2G  Ae^j  6a 

Ao  =  3K  Ae  -  3K  Ac ' '  6b 

If  we  define  the  second  term  of  the  right  hand  side  of  Equation  6  as  an 
equivalent  inelastic  deviatoric  and  volumetric  stress  increment,  we  then  can 
rearrange  Equation  6  in  the  following  form: 

ASy  +  ASy  =  2G  Aey  7a 

Ao  +  Ao’  1  =  3K  At'.  7b 

Where: 

i» 

ASjj  =  deviatoric  stress  increment 

Ao’’  =  volumetric  inelastic  stress  increment 

The  incremental  stress-strain  relation  in  terms  of  the  element  coordinates 
can  then  he  obtained  by  adding  Equations  7a  and  7b.  Thus 

Ao^  +  (AS-jJ+SijAo")  =  2G  Aeij  +  3K  (Sij  Ae  8 


Where: 
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An,-  j  I'lns!  if*  :>•  i i ':>!!  inf. ii'ft'irrvl  In  the  oHumiil  onnnHn.ihvi 
In  matrix  form,  tlio  incromentol  stress  strain  relation  is  then  given  by 
the  following  relation,  whore  the  axis  directions  are  dcfincxl  In  Figure  l. 
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9a 


Where: 


K  +  !yr~  G 


9b 


D2  =  K  -  £  G 


D3  =  2G 


9c 


9d 


For  plane  stress  conditions,  the  following  boundary  conditions  are  known: 
Ao22  =  Ao13  =  Aa23  =  Ai;12  =  Ac23  =  0  1C 


Hence  the  incremental  stress-strain  relation  is  given  by: 
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Aon 
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Ao33 
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to33 
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A°i3 
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Di^/Di 
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Where: 
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and  the  strain  component  in  the  normal  direction  is  given  by: 


Ae 


22 


=  Ao22  ~  P2  (Avyti  +  A(.33) 


11a 
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Ao22 
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D1 
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IV  m  1 1 1 .  i  i  n  .'it  1 1 1  j  1 1  <  t  n  m  I  i  I  i<  m:;#  tin*  l-i|  lowing  I  x  min  1, 1 1  y  ntinlil  inn::  .m*  kmr.,n: 


A°12  =  A°23  ~  A'  22  “  A'  12  "  A'  23  ~  0 

Ilcncc,  the  fo.Uowi.mj  incremental  stress-strain  relation  results: 


The  normal  stress  in  the  third  direction  can  be  computed  from: 

I  I 

An22  "  ^A'  1 L  '  A'  ~  A°22  ^ 

Equations  11,  13,  and  14  have  been  implemented  in  the  computer  code  to 

calculate  the  stress  increments  for  plane  strain  or  plane  stress  conditions. 

Por  a  prescribed  vector  of  incremental  strains,  the  corresponding  elastic  incre¬ 
mental  stress  vector  can  Ixv  computed  from  these  equations  once  the  elastic 
stiffness  coefficients  and  the  inelastic  stress  increment  vector  have  been  cal¬ 
culated  from  the  endochronic  parameters. 


2.2.2  Linearization  of  Endochronic  Formulation 

The  incremental  stress-strain  relation  given  by  Equation  9a  is  expressed 
in  terms  of  an  elastic  and  an  inelastic  stress  vector.  This  relation  is 


adequate  when  the  concrete  stiffness  matrix  does  not  have  to  be  combined  with 
the  crack  constitutive  relation.  For  this  cases,  the  endochronic  stress 
strain  relation  needs  to  be  formulated  in  the  following  form: 


{ao}=  [d']  {Ac} 


Where: 


|ao|  =  incremental  stress  vector  referred  to  the  element  coordinates. 
|A'e^  =  prescribed  incremental  strain  vector  referred  to  the  element 


coordinates . 


D  =  matrix  of  elastic  stiffness  coefficient  for  linearized  endochro¬ 


nic  formulation. 
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In  the  computer  program,  the  above  calculations  are  performed  only  -or  the 
in  plane  normal  and  tangential  stresses  and  strains,  and  for  the  normal  stress 
and  strain  in  the  direction  perpendicular  to  the  plane  considered.  Thus,  once 
che  deviatoric  stress  and  strain  comfx>nents  for  each  direction  is  calculated 
the  coeficients  j  are  computed  together  with  its  volumetric  component. 

bet  the  variable  x  be  defined  by: 

X,H  =  F  +  3K  (£•  L  +  V  .  l') 

3  zIf 

Then,  the  elements  in  matrix  [d' 3  are  given  by  the  following  equations 
for  the  general  case  of  plane  stress  or  strain: 

d'  (1,1)  »  K  +  i  G  -  .°iI.  -  Xn  /Bii  -  Bnn> 

3  3Z2h  (  T~ 
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d'  (2,2)  =  K  +  i  G  -  -  X22  ,B22  ~  V 

3  3Zyh  (  3  ' 

D*  (3,3)  »  K  +  1  G  -  oil  -  X33  B33  -  Bnn. 

3  3Z2h  (  ~3~) 

D  (4,4)  =  2G  -  X44  B44 

D*  (1,2)  =  K-  2g_oH  _x  b  _  Bnn 
3  3?2h  11  (  22  — ) 

d’  (1,3)  =  K-2  G  -  £H 


3S2h  Xu  <%3 '  ¥> 


D  (1,4)  =  -XU 

d'  (2,1)  =  K  -  2G  -  nil  -  X22  (BU  -  V 


3Z2h 


~T‘ 

ni 

3 


d'  (2,3)  =  K  -  -  i’ll  -  X22  (B-j _  -  Bnn 

3  3?2h  2  33  ““ 

D  (2,4)  “  “  X22  B^4 


D  (3,1)  =  K  -  2G  -  ol!  -  B, ,  -  Bn„ 

3-  3Z2h  33  <  1  -K) 

*  3 


23a 

23b 

23c 

23d 

23e 

23f 

23g 

23h 

23i 

23  j 
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d’(3,4)  =  -X33  D44  23m 

D*  (4,1)  =  -x44  {n11-  23n 

U* (4,2)  =  -  X44 j  n22  -  nnn  )  23o 

D*  (4,3)  =  -  X44  (B33  -  )  23p 

Where: 

Bnn  *  D11  *  °22  +  n33  2M 

The  stress  strain  relation  given  in  Equation  15  is  referred  to  the 

element  local  coordinate  system,  which  for  the  computer  program  has  been  assumed 

to  be  oriented  along  the  principal  stress  axis.  To  obtain  the  incremental  stress 

strain  relation  in  terms  of  global  coordinates  the  stiffness  matrix  [D  3  shall 

be  transformed  to  the  global  axis  my  means  of  the  following  relation: 


[»']„  =[t]T  W]  [t]  : 


Where: 


[ D  ]  g  =  stiffness  matrix  for  concrete  in  global  coordinates 

[t  j  =  transformation  matrix  given  by: 

r  _ 2  r2  1 


H  * 


(■•2-02 


C  =  cos  it 


S  =  sin  a 


a  =  angle  between  global  and  local  coordinate  axis 
(See  Figure  1) 


2.3  Constitutive  Relations  for  Reinforcement: 

The  constitutive  relations  for  the  reinforcement  subjected  to  monotonic 
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reinforced  concrete  structures.  A  typical  stress-strain  curve  for  the  reinforce- 
nent  subjected  to  monotonic  loading  is  shewn  in  Figure  2.  Three  different 
stages  cf  behavior  are  evident,  namely,  the  elastic  range,  the  plastic  range, 
and  the  strain  hardening  range.  During  tho  elastic  stage,  the  relation  between 
stress  and  strain  is  linear  and  is  given  by  the  modulus  of  elasticity  of  the 
reinforcement.  For  the  plastic  range,  the  strain  increases  continously  at  a 
constant  stress  and  the  modulus  of  elasticity  is  zero.  For  the  strain  harde¬ 
ning  region,  a  nonlinear  relation  exists  between  stress  and  strain,  and  a  much 
more  complicated  stiffness  relation  has  to  be  determined  from  experimental 
data. 

The  following  relations  have  been  suggested  (7  )  to  model  the  stress- 
strain  curve  for  monotonic  loading  up  to  failure. 


Elastic  Region  (t'<t-y)  j 

fs  -  E‘- 

26a 

F,s  =  29000  Ksi 

26b 

Plastic  region  (e y<« 

sh) 2  fs  =  fy 

26c 

E  =  0 

2Gd 

Strain  hardening  region  (i  sh<!  <|  sU) s 

Where: 

fg  =  steel  stress  at  strain  >. 
fy  =  steel  yield  stress 
fgu  =  steel  ultimate  stress 
t  =  actual  steel  strain 
r  =  Steel  yield  strain 

(  =  steel  strain  at  the  initiation  of  strain  hardening 
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mi 

F.  =  Elastic  mcxlulus  of  elasticity 
s  ' 

E  =  Modulus  of  elasticity  at  given  strain 

The  nonotonic  stress-strain  curve  serves  as  an  envelope  for  specimens 
subjected  to  rojxvitod  loading.  li| x >n  initial  loading,  the  stress-strain  is  si¬ 
milar  to  that  for  nonotonic  loading.  Upon  unloading,  the  stiffness  is  similar 
to  the  linear  loading  stiffness  but  a  residual  displacement  will  be  observed  if 
the  specimen  has  been  strained  to  the  plastic  range.  When  the  specimen  is  sub¬ 
sequently  loaded,  the  st.-iss-s  train  relation  is  linear  until  it  coincides  with 
the  nonotonic  stress-strain  curve,  whereupon  it  follows  the  virgin  stress- 
strain  relationship. 

The  above  constitutive  relations  are  valid  for  uniaxial  stress  states 
only.  For  reinforcement  oriented  along  3  arbitrary  directions,  the  constitutive 
relation  is  given  by: 


Hs  =  [d]s  m 


Where: 


3 

•A o  ■  =  incremental  stress  in  reinforcement  along  bar  orientations 

. 

■  ' 

.  Af  .  =  prescribed  incremental  strain  in  reinforcement  along  bar 

.  J  * 


orientations 


[d]S  =  reinforcement  stiffness  matrix 


The  reinforcement  stiffness  matrix  referred  to  the  bar  directions  as  given 


Where: 


P^,P2»Ps  =  reinforcement  ratios  along  bar  directions  1,  2,  and  3. 

E  E-,E  =  modulus  of  elasticity  of  reinforcement  along  bar  directions 
1*3 

1,  2,  and  3. 
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orientation  used  to  compute  the  a  mere  (  o  stiffness  matrix.  In  this  case ,  tlie 
stiffness  nvatrix  for  tlx.*  reinforcement  must  be  transformed  to  t'no  principal 
axis  by  the  following  relation: 


hs  -  r 
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Wlie  re: 


r  Is 

J^DJ  =  stiffness  matrix  for  the  reinforcement  referred  to  the  principal 
axis  . 

•Jins  incremental  stresses  in  the  reinforcement  can  lx:  calculated  once  tlie 
prescribed  incremental  strains  in  the  reinforcement  and  the  modulus  of  elasticity 
for  each  bar  direction  has  been  determined.  The  modulus  of  elasticity  for  each 
bar  direction  required  by  Equation  28  is  obtained  from  the  stress-strain  curve 
for  the  reinforcement  given  in  Equation  26.  It  must  be  noted  however,  that 
when  the  strain  in  the  reinforcement  is  larger  than  tlie  yield  strain  t  and 
the  element  is  unloaded,  a  new  yield  point  must  be  defined  at  the  maximum  strain 
achieved  during  the  loading  step.  The  computer  subroutine  that  calculates 
the  incremental  steel  stresses  determines  whether  the  bars  are  being  loaded  or 
unloaded  and  computes  tlie  stiffness  fc  each,  bar  direction  according  to  the 
prescribed  total  strains.  For  bar  directions  that  are  unloaded  the  code  auto¬ 
matically  shifts  the  position  of  the  yield  strain  to  obtain  the  correct  stiff¬ 
ness  for  the  reinforcement  direction  considered. 

2 . 4  Constitutive  Relations  for  Postcracking  Shear  Transfer  Mechanisms: 

The  postcracking  shear  transfer  mechanisms  ha\«  been  identified  as  tlie 
interface  shear  transfer  on  the  rough  surfaces  of  the  crack  and  the  dcwel  action 
of  the  reinforcement  crossing  the  crack.  The  interface  shear  transfer 
mechanism  is  used  to  describe  both  the  bearing  and  frictional  forces  generated 

at  open  cracks  as  the  protuding  particles  on  each  side  of  the  cracked  surface 
ccme  into  contract.  Various  experimental  investigations  (  9  ,11  ,15  )  have 
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shear  transfer  stiffness  are  tlie  initial  crack  width,  the  axial  sKffness  of 
the  reinforcement  crossing  tlic  crack,  and  the  application  of  cyclic 
loading. 

The  dcwel  action  mechanism  is  provided  primarily  by  the  herding  and 
shearing  stiffness  of  the  reinforcement  as  a  tangential  displacement  is  experi¬ 
enced  along  the  crack  length.  The  dowel  stiffness  of  the  bar  depends  mainly 
on  the  bar  diameter,  the  concrete  tensile  strength,  the  axial  stress  in  the 
reinforcement  and  tlx.'  application  of  cyclic  loadings. 

On  an  cracked  surface  of  a  reinforced  concrete  element,  both  mechanisms 
are  activated  simultaneously  to  transfer  the  applied  shear  force  across  the 
crack.  A  complete  mathematical  description  of  the  forces  and  displacements 
cxiorienced  across  the  crack  can  bo  obtained  if  a  flexibility  relationship  of 
the  following  form  can  be  established: 


Where : 

=  normal  displacement  at  crack 

AS  -  tangential  displacement  at  crack 

Aon  =  normal  stress  at  crack 

Aont  =  tangential  stress  at  crack 

Fl'  F2,  F3,  F4  =  flexibility  coefficients  at  crack 

A  =  area  of  shear  plane 
c 

The  flexibility  coefficients  required  for  Equation  30  have  been  derived 
in  Ref.  12  by  applying  incremental  unit  normal  and  shearing  stress  as  at  the 
cracked  surface  and  calculating  the  associatted  incremental  normal  and  shear 
displacements . 

Coefficient  F^  reflects  the  change  in  normal  displacement  experienced  at 


If. 
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This  coefficient  can  be  simply  described  by  the  inverse  of  tl>e  normal  restraint 
stiffness  provided  by  the  reinforcement  crossing  the  crack,  defined  by: 

Ik 


Hence,  coefficient  F^  is  given  by: 


32 


n  8240db4 

Where  the  normal  restraint  stiffness  1^,  is  calculated  from  the  relation  pro¬ 
posed  by  Jim6ncz,  et  al.  (13) 

The  flexibility  coefficient  F2  represents  the  increase  in  normal  displace¬ 
ment  caused  by  the  applied  shear.  If  it  is  assumed  that  the  increase  in  crack 
width  or  normal  displacement  is  caused  mainly  by  the  interface  shear  transfer 
stresses,  ttie  increase  in  crack  width  can  be  calculated  from  the  normal  stresses 
induced  by  the  applied,  shear.  Based  on  an  expression  proposed  by  Reinhardt 
and  Walraven  ( 17  ) ,  the  change  in  crack  width  can  be  calculated  from: 


F2  = 


Where: 


a 


-0.63  -.552  ,  •  1 

0.176c  +  (0.22c  -  1.034)  f„ 


0.1353c  "°*8  +  (0. 164c” *?07-l. 379)  f'c 
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=  ratio  of  interface  shear  transfer  stiffness  to  the  sum  of  the 
interface  shear  transfer  and  dowel  action  stiffness. 

I 

f  =  concrete  compressive  strength  (ksi) . 
c 

c  -  initial  crack  width  (in.). 

The  flexibility  coefficient  can  be  calculated  if  it  is  assumed  that  the 
increase  in  shear  displacement  is  caused  by  the  reduction  in  the  interface 
shear  transfer  stiffness  associated  with  a  larger  crack  width.  In  mathematical 
terms. 
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tlxj  change  in  crack  width  will  be  obtained  once  the  equation  for  tho  cocf ficiont  F^ 
is  obtained.  Note  that  in  Equation  34,  the  change  in  normal  displacement  with 
respect  to  the  change  in  normal  stress  is  proportional  to  the  normal  restraint 
stiffness.  Thus,  liquation  34  can  bo  rewritten  as: 

f3  -  if1-  •  35 

Kn  dK^c 

The  flexibility  coefficient  represents  the  incremental  shear  displacement 
experienced  at  the  crack  when  an  incremental  unit  tangential  shear  force  is 
transferred  across  the  crack  The  shear  displacement  is  inversely  proportional 
to  the  total  stiffness  provided  by  the  interface  shear  transfer  and  the  dcwel 
action  mechanisms.  Given  the  stiffness  of  both  mechanisms,  the  function  can 
be  calculated  from  the  following  equation: 


+  *d 


Where: 


Ka  =  interface  shear  transfer  stiffness  (K/'in) 

=  dcwel  stiffness  of  reinforcement  crossing  the  crack  (k/in) 

Based  on  a  review  of  several  relations  available  for  the  interface  shear 
transfer  and  dowel  action  stiffness,  the  following  equation  was  selected 
from  Reference  12. 


Interface  Shear  Transfer  Stiffness: 


Ka  " 


3.9(c-0.002)  +  1.09x10”'  (3.4  x  10  -  K 


Where: 


K_  =  interface  shear  transfer  stiffness 

a 

=  normal  restraint  stiffness 
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n2  * 

-54  for  Vd 

>  °‘9vdu 

dowel  displacement  (in) 
axial  stress  in  reinforcement  (Ksi) 
yield  stress  of  reinforcement  (ksi) 
ultimate  dowel  capacity  of  reinforcement  (K) 
bar  diameter  (in) 

The  ultimate  dowel  capacity  of  the  reinforcement  is  controlled  by  whether 
the  dowel  will  fail  by  yielding  of  the  reinforcement  or  by  concrete  splitting. 
Thus,  the  ultimate  dowel  load  is  given  by  the  smaller  of  the  values  predicted 
by  the  following  relations: 

Failure  by  yielding  of  the  reinforcement: 

2 


Vdy  *  0.92^  J~t ~  t'c 
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Where: 


=  dowel  failure  load  caused  by  yielding  of  the  reinforcement  (Kips) 


j 

5 


f  -  ylt '  1  <1  nr  Ihn  tvinlm •  mi'iil  (Kni) 

y 


Failure  by  concrete  splitting: 
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Where: 

V,  «  dowel  failure  load  caused  by  concrete  splitting  (Kips) 
do 

b  *  net  width  of  section  perpendicular  to  load  direction  (in) . 
n 

■  number  of  bars  per  layer. 

c  «  smaller  of  side  or  bottom  concrete  cover  of  the  reinforcement  (in) 
m 

Thus,  the  flexibility  coefficient  can  be  obtained  fran  Equation  36 

once  the  interface  shear  transfer  and  dcwel  action  stiffness  have  been  calcu- 

1  .* 

lated  from  Equations  37  and  38. 

The  equation  for  coefficient  F3  can  new  be  presented  once  the  first 
derivative  of  Equation  37  with  respect  to  the  initial  crack  width  is  computed: 


F3  - 
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Where: 


av  «  shear  stress  increment  applied  in  previous  step, 
o 
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2,5  Constitutive  Relations  for  Cracked  Reinforced  Concrete 

The  incremental  stress  vector  induced  by  a  prescribed  strain  increment  in 
a  reinforced  concrete  element  can  be  obtained  from  the  incremental  stress 
vectors  sustained  separately  by  the  concrete  and  the  reinforcement,  provided 
that  the  concrete  element  has  not  cracked.  If  we  assume  that  the  average 
strains  in  the  concrete  and  the  reinforcement  are  equal,  then  the  total  incre¬ 
mental  stress  can  be  calculated  from: 

(to}T  -  [fn]  ♦  [»]*]  {4,} 


46 


Whoio: 


|Aa|r  =  total  incremental  stress  vector 

jAf.j  =  prescribed  incremental  strain  vector 

J^D*  =  stiffness  matrix  of  uncracked  concrete  element 

r  i  s 

I)  -  stiffness  matrix  of  reinforcement 
For  reinforced  concrete  elements  where  the  principal  tensile  stress  has 
exceeded  the  maximum  tensile  strength  of  the  concrete,  the  incremental  stress 
vector  is  a  function  of  the  tangential  and  normal  stresses  transferred  across 
the  crack.  For  this  cases  the  constitutive  relation  given  in  Equation  46  has 
to  be  modified  as  described  subsequently. 

When  the  principal  tensile  stress  exceeds  the  maximum  tensile  strength  of 
the  concrete  the  prescribed  incremental  strain  calculated  for  the  current  step 
has  to  be  divided  into  the  incremental  strain  required  for  the  element  to  crack 
and  the  remaining  incremental  strain  necessary  to  complete  the  total  incremental 

i 

strain  computed  for  the  current  step.  Hence, 

{4=  {is}1*  {ic}2  47 

Where: 

|Aej  =  total  incremental  strain  at  current  time  step 
|  Ae| ^ =  incremental  strain  required  for  crack  initiation 
|  Aej  2  =  incremental  strain  required  to  complete  the  total  strain  increment 
assigned  to  current  step. 

The  incremental  strain  required  for  the  element  to  crack  is  estimated  from 
the  proportion  of  the  stress  increment  at  which  the  principal  stress  equals  the 
tensile  strength  of  the  concrete.  Said  proportion  is  given  by: 


p  *  ft  -  si 


Where: 
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=  Proportion  of  stress  increment  at  which  crac'cing  ocurred 
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tensile  St  length  of  (•onri  el  n 
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n-1  = 


=  principal  tensile  stress  nt  current  step. 

principal  tensile  stress  at  previous  step 

Thus,  the  incremental  strain  at  which  cracking  ocurred  is  given  by: 

|  Ac}1  =  P{A'} 

The  remaining  incremental  strain  to  be  applied  to  the  cracked  element 
during  the  current  step  is  then  given  by:  ' 

| Ac}2  =  (  1  -  P)  {A.  } 

The  incremental  strain  applied  to  the  cracked  element  required  to  complete 
the  current  step  is  distributed  between  the  cracked  and  uncracked  sections  of 
the  element  according  to  the  following  equations: 

|  At  }  cr  =  a  |  At :}  2 

(1-a  )  { Ae} 2 

strains  contributed  by  the  cracks  in  the  element 
strain  contributed  by  the  uncracked  section  of  concrete  within 

i 

the  cracKS . 

a  =  proportion  of  incremental  strain  provided  by  cracks  within  the 
element 

Onoe  the  strains  contributed  by  the  cracks  are  known  the  average  normal 
and  tangential  displacements  can  be  determined  from  the  crack  spacing,  as 
given  by  the  following  relations: 
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Where: 
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incremental  strain  normal  to  crack  contributed  by  the  cracks 
incremental  strain  parallel  to  crack  contributed  by  the  cracks 
crack  spacing 
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II  should  lx>  noli'il  th.il  in  Ii'Iiim;  oi  llu*  lour  si  ross  oonii'out'nl  n  t'onsitli'rod 
in  I  ho  computer  c'olo,  llu1  lonsl  il  ul  ivu  ndntiou  for  the  cracks  pn\sontoil  in 
Initiation  30  con  new  )m  expressed  in  to  nan  of  the  crack  strains  by  the  fol  lowing 
ro la l ion: 
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The  strains  contributed  by  the  uncracked  concrete  between  the  cracks 
are  used  to  determine  the  endochronic  parameters  required  to  calculate  the  con¬ 
crete  stiffness  matrix  given  in  Equation  15.  As  the  stresses  in  the  uncracked 
concrete  and  in  the  crack  have  to  be  similar,  the  proportion  of  incremental 
strain  taken  by  the  uncracked  concrete  (See  Equation  51)  is  determined  in  an 
iterative  fashion  from  the  following  relation  once  the  crack  flexibility  matrix 
and  the  concrete  stiffness  matrix  are  known: 


If  the  new  value  for  the  uncracked  strains  are  within  tolerable  limits  of 
the  assumed  uncracked  strains,  then  the  convergence  requirement  that  the 
stresses  in  the  solid  concrete  be  equal  to  the  stresses  in  the  cracked  concrete 
has  been  satisfied.  Otherwise  the  previous  value  of  uncracked  strains  is 
replaced  by  the  latest  vector  of  uncracked  strains,  a  new  cracked  strain  vector 
is  computed  and  the  cracked  flexibility  matrix  and  the  uncracked  concrete 
stiffness  matrix  are  calculated  again.  The  iterative  procedure  is  continued 
until  the  convergence  requirement  is  satisfied. 

The  incremental  stress  for  the  cracked  elements  can  then  be  calculated 
from  the  following  constitutive  relation  for  cracked  reinforced  concrete. 


;n 
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Thus,  the  total  incremental  stresses  in  the  concrete  attained  during  the 
current  step  are  calculated  from: 

=  {An}1  +  {am}^  55 

=  total  incremental  stresses  in  concrete  during  current  step 
=  incremental  stresses  required  to  crack  the  element 
=  incremental  stresses  for  the  cracked  element 
The  total  stresses  sustained  by  the  reinforced  concrete  element  is  obtained 
by  adding  the  steel  stresses  to  the  concrete  stresses  computed  from  Equation  55. 

It  should  be  noted  that  the  crack  formation  criteria  used  in  the  computer 
code  is  based  on  the  maximum  tensile  stress  criteria.  A  tensile  crack  is  formed 
whenever  the  principal  tensile  stress  exceeds  the  maximum  tensile  strength  of 
the  concrete.  Once  the  crack  Is  fomed,  the  stress  in  the  concrete  normal  to 
the  crack  is  set  to  zero  and  the  concrete  strain  at  which  the  crack  ocurred  is 
stored.  The  crack  is  assumed  to  close  whenever  the  concrete  strain  is  smaller 
than  the  strain  at  which  the  crack  ocurred.  For  closing  cracks,  the  constitutive 
relations  used  are  similar  to  those  used  for  the  initially  uncracked  concrete. 

The  crack  is  assumed  to  open  again  whenever  the  concrete  compressive  stress 
drops  to  zero, whereupon  for  subsequent  loading,  the  constitutive  equations  for 
cracked  reinforced  concrete,  given  by  Equation  54,  are  used. 


3.  Computer  Program  for  Nonlinear  Analysis  of  Reinforced  Concrete 
3 . 1  Introduction: 

The  computer  program  developed  to  model  the  nonlinear  behavior  cf  reinforced 


concrete  plane  stress  or  plain  strain  elements  based  on  the  theoretical  concepts 

established  in  Chapter  is  given  in  Appendix  A  2.  A  flew  chart  is 
first  discussed  herein  to  establish  the  sequence  of  principal  oi«rations 


rFnTTri'"'!  •’S'Wpr  •wr-r»p'r™ . . 
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l*'rlonu!»l  l>y  tin*  piixii.iin,  lolhMin  by  .1  doner i pi  ion  ol  the  .lrt.ivll.lon  imrlnnrMl 
within  each  mil  until  inn  included  in  the  «,odo. 

'Ilio  otmputor  Bubrnufinoft  prortonl  <xl  In  A|H*’ndix  A2  arc  lo  Im  inmnxirntcd 
in  an  existing  conjxitor  prrxir.nn  nt  tlx?  Air  Force  Whalens  laboratory  identified 
by  the  acronym  of  RAMSON  (6) .  This  cock;  is  used  to  perform  nonlinear  dynamic 
it.ialys.is  of  plane  and  axisynmotric  problems  but  at.  present  considers  only  the 
nonlinear  stress-strain  relation  for  concrete.  The  material,  subroutine 
discussed  subsequently  should  enhance  the  analytical  capabilities  of  the 
computer  code  SAMSOM. 


Hie  main  purpose  of  the  computer  code  is  to  calculate  the  incremental 
stresses  induced  by  a  prescribed  vector  of  incremental  strains  for  an  uncracked 
or  cracked  reinforced  concrete  element.  The  material  subroutines  require  that 
the  prescribed  incremental  strain  vector  be  defined  beforehand.  The  code  does 
not  include  an  equation  solving  subroutine  as  those  operations  will  be 
performed  by  tne  principal  code  SAMSOM. 
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(’ln'ok  IT  ik<w  iT.*+.fi  have  h  >n»xl.  If  lltt*  priivijMl  tensile  Ml  less 
is  smnl  lor  than  the  concrete  tensile  strength  go  to  stop  12. 

5.  Determine  order  and  numlx.'r  of  cracks. 

6.  Update  stresses  in  current  stop  to  state  of  incipient  cracking. 

7.  Start  iteration  for  uncrackcxl  strains.  Determine  distriljution 
of  incremental  strain  vector  loft  in  current  step  between  the 
uncracked  and  cracked  concrete. 

8.  Determine  total  cracked  directions. 

9.  Determine  crack  flexibility  matrix. 

10.  Determine  uncracked  concrete  stiffness  matrix. 

11.  Compute  new  vector  of  uncracked  strains. 

12.  If  number  of  iterations  for  uncracked  strains  is  smaller  than 
three  go  to  step  D7. 

13.  Compute  stresses  in  cracked  concrete  element. 

E.  Compute  incremental  stress  vector  for  reinforcement  caused  by  pres- 
scribed  incremental  strains. 

F.  Transform  incremental  stress  vectors  for  concrete  and  reinforcement 
to  global  directions. 

G.  Update  stresses  and  strains. 

H.  Proceed  to  next  strain  increment 

I.  End 

3 . 3  Description  of  Code  Subroutines : 

The  computer  code  presented  in  A^jondix  A2  contains  14  subroutines  in 
addition  to  the  main  section  of  the  program.  The  main  section  is  used  to 
compute-  the  incremental  strain  vector  according  to  the  analysis  desired 
and  to  read  the  control  and  material  data  required.  The  subroutines  are 
described  subsequently. 


Subroutine  (Y)KF:  This  nubmiil  ine*  <m1  led  liy  IIm'  main  ion  of  Mw 
program,  computes  tlxj  constant  coefficients  required  for  the  endtx:hrnn.vo  ctf na¬ 
tions  given  in  Apjxsndix  Al.  Input  requ:  rod  are  the  concrete  canpressivo 
strength,  the  reinforcement  yield  stress  and  the  reinforcement  ratios  in  each 
direction.  This  subroutine  iVxs  not  call  n.:y  other  subroutine. 

Subroutine  CRACUKt  This  subroutine  is  called  by  subroutine  MATERI  to 
determine  if  initial  cracks  have  ocurred  during  the  current  stress  increment  or 
to  determine  if  the  cracks  in  a  previously  cracked  element  are  closed  or  open. 
Input  required  are  the  stress  and  strain  vcctirs  for  the  current  and  past  steps 
and  the  strain  at  which  the  crack  opened  previously.  This  subroutine  does  not 
call  any  other  subroutine. 

Subroutine  CRASH:  This  subroutine  is  called  by  subroutine  ONECRA  to 
compute  the  crack  flexibility  matrix.  Inquired  input  are  the  previous  strain 
and  stress  vector,  the  proportion  of  incremental  strain  contributed  by  the 
cracks,  and  the  geometric  proper  ties  of  the  element  such  as  bar  diameter, 
concrete  cover,  number  of  bars,  etc.  This  subroutine  does  not  call  any  other 
subroutine. 

Subroutine  FUNEND:  This  subroutine  is  called  by  subroutines  MATERI 
and  ONECRA  to  compute  the  endochronic  model.  Subroutine  MATERI  calls  FUNEND 
to  compute  the  incremental  stress  vector  for  uncracked  concrete  caused  by  the 
prescribed  incremental  strains.  Subroutine  ONECRA  calls  FUNEND  to  compute  the 
uncracked  concrete  stiffness  matrix  for  the  uncracked  incremental  strain 
vector  and  to  compute  the  cracked  concrete  incremental  stresses.  Required 
input  is  the  stress  and  strain  vector  for  the  current  and  previous  step  and 
whether  the  element  is  in  plane  strain,  plane  stress  or  axisynmetric.  This 
subroutine  calls  subroutines  INVAR  and  PRIN. 

Subroutine  GLOBAL;  This  subroutine  is  called  by  subroutine  MATERI  to 
transform  the  conputad  concrete  and  reinforcement  stress  vectors  in  the 


21 


current  step,  rofonrxl  to  Llx?  principal  axis,  Kick  to  the  qlnlvp  mordi  tinier?. 
Required  input  is  the  stress  vector  to  lx?  computed  and  tlie  angle  between  tiie 
principal  axis  and  the  global  coordinate  axis.  This  subroutine  docs  not  call 
any  other  subroutine. 

Subroutine  INV:  This  subroutine  is  called  by  other  subroutines  to  compute 
the  inverse  of  a  given  matrix.  Required  input  are  the  matrix  to  be  inverted, 
the  array  where  the  results  will  be  stored,  and  the  order  of  the  matrix.  This 
subroutine  does  not  call  any  other  subroutine. 

Subroutine  I  WAR:  'Hits  subroutine  is  called  by  subroutine  FUNEND  to 
compute  the  stress  and  strain  invariants  for  the  current  values  of  the  strain 
and  stress  vectors.  Input  required  are  the  current  values  of  the  strain  and 
stress  vectors,  together  with  the  incremental  strain  vector.  Subroutine 
INVAR  does  not  call  any  other  subroutine. 

Subroutine  MATER1:  This  subroutine  is  called  by  the  main  section  of  the 
program  to  compute  the  incremental  stress  vector  in  the  concrete  and  reinforce¬ 
ment  caused  by  the  prescribed  strain  vector.  Input  required  for  this 
subroutine  is  the  current  and  previous  step  vectors  of  stress  and  strain,  the 
total  stress  vector  of  steel  stresses  and  whether  the  element  is  in  plane 
stress,  plain  strain  or  axisymnetric.  This  subroutine  calls  the  subroutines 
ROTATE,  FUNEND,  CRACHK,  ONECRA,  and  STEEL. 

Subroutine  MA*IMU;  This  subroutine  is  called  by  other  subroutines  as 
required  to  multiply  two  given  matrices.  Input  required  is  the  name  of  the  two 
matrices  to  be  multiplied,  the  name  of  the  array  where  the  results  will  be 
stored,  and  the  order  of  the  matrices.  This  subroutine  does  not  call  any  other 
subroutine. 

Subroutine  ONECRA:  This  subroutine  is  called  by  subroutine  MATERl  when 
initial  cracking  in  an  element  is  detected.  The  subroutine  updates  the  stresses 
to  the  state  of  incipient  cracking,  determines  the  preportion  of  total  strain 
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iteration  on  the  cracked  oik!  uneraekcd  concrete  strains,  and  confutes  the 
stresses  in  the  cracked  concrete.  'Hie  input  required  are  Die  current  and 
previous  step  strain  and  stress  vectors  referred  to  principal  and  global  di¬ 
rections,  the  strains  at.  which  cracking  previously  ocurred,  tire  angle  between 
the  principal  and  global  directions  and  the  cracking  direction.  This 
subroutine  calls  subroutines  CRASTI  and  RINEND. 

Subroutine  PRIN:  This  subroutine  is  called  by  subroutines  ROTATE  and 
FUNEND.  Subroutine  ROTATE  calls  PRIN  to  determine  the  principal  stresses  or 
strains  from  the  global  stress  or  strain  vector.  Subroutine  FUNEND  calls 
PRIN  to  arrange  tl»e  maximum,  intermediate  and  minimum  principal  stress  in 
increasing  order  of  magnitude.  Input  required  is  the  stress  or  strain  vector 
and  whether  principal  stresses  are  to  be  computed  or  not.  This  subroutine 
does  not  call  any  other  subroutine. 

Subroutine  ROTATE:  This  subroutine,  called  by  MA.TER1,  carputes  the  prin- 
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cipal  stresses  and  strains  for  the  current  and  previous  strain  and  stress  vec¬ 
tors.  Input  required  are  the  strain  and  stress  vectors  for  the  current  and 
past  steps.  This  subroutine  calls  subroutine  PRIN. 

Subroutine  Steel:  This  subroutine  is  called  by  MATERl  to  compute  the 
incremental  stress  vector  induced  in  the  reinforcement  by  a  prescribed  strain 
increment.  Required  input  are  the  current  and  previous  strain  vectors  in 
principal  and  global  coordinates .  It  should  be  noted  that  in  this  subroutine 
the  bars  are  assumed  to  he  oriented  along  the  element  global  directions  (See 
Figure  1),  This  subroutine  does  not  call  any  other  subroutine. 

Subroutine  ZER:  This  subroutine  is  called  by  the  other  subroutines  as 
required  to  initialize  the  values  of  arrays  to  zero.  Input  required  is  the 
matrix  to  be  initialized  and  its  order.  No  other  subroutine  is  called  by 


ZER. 


4.  Ntmi'tio.i!  UtmiiIIm: 

4.1  Inlmduel  Jon: 

In  Mils  (iuipUn  ,  Mm>  u'nuIMs  ol  ,« ;< *v* * i . i I  px|x‘i  imrnl ul  pioqiuiMM  .no 
with  Uic  results  pi.vtU.oUxl  by  the  ixnipulm:  prixjrani  in  older  to  ascertain  its 
applicability.  Connin’. sons  arc  made  with  exiieri  mental  stress-strain  curves  for 
plain  concrete  in  uniaxial  and  biaxial  stress  states,  and  for  reinforced  concrete 
panels  subjected  to  biaxial  stress  states. 

4.2  Stress  strain  Curve  for  Plain  Concrete; 

In  this  section,  the  stress  strain  curve  for  plain  ooncrete  predicted  by 
the  computer  code  is  compared  with  experimental  results  for  uniaxial  and 
biaxial  stress  states.  The  prescribed  strains  are  increased  in  0.0001  incre¬ 
ments  and  the  resulting  induced  stresses  are  computed  from  the  constitutive 
relation  given  in  Equation  15  for  the  corresponding  stress  states. 

In  Figure  3,  the  computed  stress  strain  curve  for  uniaxial  loading  of 
plain  concrete  is  compared  with  the  experimental  data  given  in  Reference  10 
for  various  concrete  strengths.  The  endochronic  formulation  coded  in  the 
program  predicts  curves  that  are  very  similar  to  the  observed  behavior. 

Figure  4  compares  the  stress  strain  curve  for  biaxial  stress  states  com¬ 
puted  by  the  program  with  the  experimental  points  determined  from  Reference  14. 

As  for  the  uniaxial  tests, the  predicted  values  agree  well  with  the  observed 
behavior. 

4  3  Stress-Strain  Curves  for  Reinforced  Concrete  Panels  Subjected  to  Biaxial 
Stress  States; 

Several  experimental  results  have  been  reported  in  the  literature  (16, 

18)  where  reinforced  concrete  panels  have  been  subjected  to  monotonic  and  cyclic 
shear  stresses.  In  this  section  the  results  of  one  test  reported  by 
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oxk'.  Additional  (xnniAiriMiim  with  oxtx'rimontal  data  should  lx?  conducted  wlicn 
tlx?  ctxlo  is  incorporated  to  tlio  main  program  to  determine  its  general  applica¬ 
bility  to  more  complex  problems. 

Tlie  test  s|xx:iiucn  tvpn-Lod  by  Fonlikaris,  et  al,  Ii.kI  a  central  section 
of  24  in  by  24  in  with  a  thickness  of  6  in  and  reinforced  in  ono  direction 
with  one  layer  of  #4  bars  at  6"  spacing  and  with  two  layers  of  #4  bars  at 
6"  spacing  in  the  perpendicular  direction.  The  concrete  compressive  strength 
was  3160  psi,  while  the  reinforcement  yield  stress  was  reported  as  61000  psi. 

The  specimen  was  subjected  to  a  monotonic  shear  stress  up  to  failure  without 
the  application  of  biaxial  normal  stresses.  The  ultimate  shear  stress 
sustained  by  the  specimen  was  reported  at  475  psi. 

The  experimental  load  deflection  curve  determined  for  this  specimen  is 
given  in  Figure  5  together  with  the  calculated  load  displacement  relation.  In 
general,  both  curves  corelate  well  but  vhe  model  fails  to  predict  the  correct 
ultimate  shear  stress.  It  has  been  observed  that  the  crack  flexibility  matrix 
becomes  jJlconditioned  at  large  strains  and  predicts  significantly  different 
results  for  the  cracked  incremental  stresses.  It  is  believed  that  the  problem 


is  caused  by  the  equation  proposed  the  flexibility  coefficient  F-. 


n , 

5.  q->nclustons: 

The  computer  code  developed  heroin  can  be  used  to  model  the  nonlinear 
behavior  of  cracked  reinforced  concrete  in  finite  element  analysis.  At  this 
point,  the  code  considers  the  nonlinear  strcss-strain  curve  of  plain  concrete 
and  of  the  reinforcement,  the  anisotropy  induced  by  complex  stress  states  ami 
by  cracking,  and  the  postcracking  shear  transfer  mechanisms.  The  code  can  be 
used  for  the  nonlinear  analysis  of  reinforced  concrete  structures  once  it  is 
incorporated  into  the  main  program  SAMSOM. 

The  program  predicts  results  that  agree  reasonably  well  with  a  limited 
amount  of  experimental  data  for  plain  and  reinforced  concrete  elements.  A 
large  number  of  additional  tests  should  be  conducted  when  the  material  subroutine 
is  incorporated  in  the  main  program  to  check  the  program  predictions  with  the 
available  experimental  data  on  reinforced  concrete  panels  reported  by 
Perdi’  aris,  et.  al.  (16)  and  by  Vecchio,  et.  al  (18).  This  effort  will  deter¬ 
mine  tue  final  application  range  of  the  material  subroutine  developed  herein. 
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b)  local  coordinates 


Figure  |:  Coordinate  Definition 


A.  1  .'•iimviry  nT  lv|U,tt  ii>iv.  foi  Kn>  It  n ‘bionic  Mi m ti *  1 

'Hks  equations  ptx'sentod  liorein  arc*  suimvirized  from  Ref.  5  ,  to  compute  the 
material  constants  and  functions  reciu.irvd  to  model  tlie  nonlinear  behavior  of 
concrete. 

Distortion  Intrinsic  'tune  Parameter  Z: 


AZ  = 


An  =  (F)  Af, 


A£  = 


F  = 


(J2  (AO  '1/2 


F  +  F 
*1  2 


Qq  (1-91) _ 

1  -  a5  [I3(0]  1/3  (1  +  g2) 


F2  = 


F2  " 


f  = 


A. 

a2  (1+ Jag  I2(o)  |  1/4  +  F5) 

=  l-a^to)  +  |a8  I2  (o)  1/4F4-a3I3(o)  [  J2  (a)  j  1/8  {1'gO 

‘  j[  +  fl  +  ^  ?3 

1  *  F2/a7 


^  = 
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J2(<  )  d+a9/n2) 
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^2  “  ^21  9  ■  J23 
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9n  =  al4  J2  1‘  /  °ned  ~  °min 

\  °max  “  a23 


h  (? 


(->med  ~  °min 
°max  ~  a23 
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Inelastic  Dil.atnncy  Parameter  X(: 
AX  =  (l)  (L)  (Af.) 


i, 

=  1-  X_ 

Xo 

A1.30 

L 

c3  r 

l“cill(0) 

‘X_  ' 

2  + 

1  21 

C2I+J2(i-) 

A1.31 

Shear  Gonpaction  Parameter  X ' : 

ax' 

”►4 

~a> 

II 

A1.32 

l' 

=  c6  (1-  X1 

/ 

Xo  ) 

A1.33 

L 

„  1/3 

°min  ^3 

1  +  |93  /  C8|3 

Al.34 

g3 

=  °'S3  |c7  “mini 

-  J 

J2(£  ) 

A1.35 

Bulk  and  Shear  Modulus: 

K 

-  1 

Eo 

A1.36 

1  +  c  x 

5 

3  (l-2v) 

G 

=  1 

Eo 

A1.37 

1  +  c5\  2 (1+v) 


Stress  and  Strain  Invariants: 

The  following  equations  are  valid  for  principal  stresses  and  .strains  only: 


X1  (o'i  =  °11  +  °22  +  J33 

A1.38 

t  In)  =  ^all~a22^2  +  <°22  "  °33)2  +  ^a33"°ll^2 

2  6 

Al-39 

I3  (a)  =  (0^)  (o22)  (^33) 

A1.40 

I^(Ae)  —  Agj_j_  +Ac22  AE33 

A1.41 

j0  (Ae)  =  (Arli."  A22)2  +  (Ac22~At33)2  +  (Ae33-Acll)2 

Al,42 

6 

=  (Ell_e22)2  +  (e22  "  e33>2  +  (^33  ”  ^11>2 

A1.43 
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r-joofJn  *  wnnoo  non 


DIMENSION  S 1U1 ( A ) » SIU 1 1 ( 1 o » 6  > » EPST  K  A  > » EPS T 1 1 ( 1 0. 6 ) 
DIMENSION  DER( 4  >  . DUM<  4  > 

COM  C  ( 60 ) » LK  1 »  A  > 

DATA  ISTRES. IESLAST. I PLAST » UGT/2 . 0 » 0 » 0 
DATA  £K1)»E2<1>»E3<1>»E<4).E5(I>/0.0»0.0»0 
DO  117  1-1.4 

EPST I1(I»J)*PREVI0U8  STRAIN  OF  ELEMENT  I  IN  DIRECTION  J 
EPSTim-CURRENT  STRAIN  IN  DIRECTION  I 
SIGIKI.  J)-PREVIOUS  STRESS  OF  ELEMENT  I  IN  DIRECTION  J 
SIGI ( I ) “CURRENT  STRESS  IN  DIRECTION  I 
TST ( I )*T0TAL  STEEL  STRESS  IN  DIRECTION  I 

EPSTI1(1»I>=0.0 
EPST I < I ) *0 
SIGIK1»I>=0 
SIGI<I>«0 
TST ( I ) =0 
.7  CONTINUE 

STINCa0«0001 

INITIALIZE  OR  READ  PREVIOUS  STRAINS  AND  STRESSES 

EPST I (4) “EPST I (4)4STINC 
IONE-1 

FC=C0NCRETE  COMPRESSIVE  STRENGTH 
FY-REINFORCEMENT  YIELD  STRESS 
ECT-MAXIMUM  TENSILE  STRAIN  IN  CONCRETE 

PS1 » PS2  »PS3  IS  REINFORCEMENT  RATIO  IN  DIRECTIONS  1.2  AND 


0131 


DATA  FC.FY.ECT.PS1.PS2.PS3/22A3, .34809. »-. 01 » .017? .0, » 
CALL  COEF(FC»FY»ECT .PS1.PS2.PS3) 


30  CALL  MATERKSIGI.SIGI1.EPSTI.EPSTI1.TST.C.DCR.DUN.ISTR 
ES) 

DO  118  1-1.4 
SIGI1(1»I)=SIGI(I> 

EPSTI1<1.I)*EPSTI<I) 

118  CONTINUE 

EPST I  <  4 ) *EPSTI <4)+STINC 
GO  TO  230 
END 

SUBROUTINE  MATER 1 (SIG.SIGl . EPS. EPS 1 »TST» C.DCR.DUN. ISTR 

ES) 

DIMENSION  ST(A) .THE (4) »DS(4.4).CRASTR(4) »DEUN(4) »DECR( 

4) 

1  =  1 

FT=C<53) 

ITE=0 
CRADIR-O 
N  TCR*0 
ITE»0 
C 

C  SAVE  STRESSES  AND  STRAINS  BEFORE  CHANGING  TO  PRINCIPAL  DIR 
ECTIONS  , 

C 

DO  119  1=1.4 
SKI  )=SIG1  ( 1 . 1 ) 

S( I )~SIG( I > 

El(I)-EPSKl.I) 

.  E< I )=EPS( I ) 

119  CONTINUE 
* 


£  CHANGE 
C 


STRESSES 


STRAINS 


PRINCIPAL  DIRECTIONS 


CALL  ROTATE(SvSl,E,El ,THE> 

I F  <  T  HE  <  2 ) , EQ . 0 , 0 . AND , THE  <  4 ) . NE , 0 , 0 ) THE ( 2 )  =  THE ( 4 ) 

DO  120  I ~ 1 »  3 

IF(CRASTR< I) .EG. 0.0)  GO  TO  120 

NCR=NCR+1 

ICR  =  I 

120  CONTINUE 

IF ( NCR .  EQ  »  0 ) GOT 02000 
IFCNCR.EQ. l.AND. ICR . EQ , 2 ) GGT02000 
IF(NCR#EQ«0)GQT02000 
GOTO  5000 

C  COMPUTE  ENDOCHRONIC  PARAMETERS  AND  INCREMENTAL  STRESSES 
2000  CALL  FUNEND <  8 »  SI , E »  E 1 »  THE , DS , I STRES » NCR ) 

C  CHECK  FOR  CRACKS  IN  CONCRETE  (MAX.  SRESS  CRITERIA)  , 
CALL  CRACHK(S»S1 »E»E1 ,CRASTR»FT »NTCR, CRADIR) 
IF(NTCR*EQ.O)GOTO  8000 

IF (NTCR* l.AND. CRADIR. EG. 2) GOTO  2400 
1F<NTCR,EQ.1. AND, CRADIR.LT, 4) GOT 02500 
IF < NTCR. EQ.2. AND. CRADIR. EQ. 4) GOTO  4000 
IF(NTCR,EQ,3)G0T0  4000 
GOTO  2500 
2400  S<2>=0 

I F ( I STRES , EQ ♦ 2 ) GOTO  8000 
GOTO  8000 

C  EXIT  IF  NEW  CRACK  IS  IN  2(THETE>  DIRECTION 

2500  WRITE <  3, 27 ) 

27  F0RMAT<5X, ' ***WARNING ! ELEMENT  HAS  CRACKED**# ' ) 

C  CRACK  IN  ONF.  DIRECTION/DETERMINE  CRACKING  DIRECTION 
ICR=CRADIR 

CALL  0NECRA(S»S1»E»E1»SIG*SIG1»EPS»EPS1»CRASTR»C»TST»T 
HE»BUN»DCR»A»FT, ICR, ISTRES) 

GOTO  8000 

8000  CALL  STEEL(ST»EPS»EPS1 »E»E1 »THE»STRMAX> 

C  CHANGE  STEEL  AND  CONCRETE  STRESSES  TO  GLOBAL  COORDINATES 
B  =  2 

CALL  GLOBAL. (ST,  A,  D) 

DO  121  1  =  1,4 
SX<I)=S<I) 

EX ( I ) =E ( I ) 

121  CONTINUE 

CALL  GLOBAL(SX»A,B> 

B-l 

CALL  GLOBAL < EX » A, D) 

DO  122  1=1,4 
EPS  < I ) =EX  < I ) 

E  < I >  =EPS ( I ) 

S< I )=SIG< I )=SX< I ) 

122  CONTINUE 
EPS  <  2 ) =E ( 2 ) 

DO  123  1=1,4 

TST (I)=TST ( I ) +ST ( I) 

SU1=SIG(I)+TST(I) 

WRITE ( 3i’8 ) SIG  < I ) » TST ( I ) , SU1 

28  FORMAT ( 5X »  ' SCON= ' , E 10 , 4 ,5X , 'STEEL® ' , E10 , 4 , 5X, ' TOSTR= ' 


,  E 1 0  »  4 ) 


FORMAT ( 5X 

CONTINUE 

STOP 

RETURN 


>/4. 

C 


R  t)  B  K  Old  INK.  I N  V  A  R  ( S  *  E  » 1 1 F  i  R  I  N 1 » S  l  N  2 » S I N  3 » D  S I N  2 » D  E I N 1 1  D  D  E! 
OPTION  base:  1 

THIS  SUBROUTINE  COMPUTES  STRESS  AND  STRAIN  INVARIANTS 
STRESS  INVARIANTS 
SAV=S( 1 )  +3 <  2 ) +S<  3 ) 

SAV=SAV/3 . 

SIN1=SAV*3, 

SIN2*-(S(l)*S(2)+S(2)*S(3)fS( 1>*S<3> ) 
SIN3=S<1)*S<2)*S<3> 

DSIN2=< <S<1)-S(2) )#*2.+<S<2)-S<3> > **2 ♦ + < S < 3 > -S < 1 > >**2. 


C  STRAIN  INVARIENTS 

DEAV=DE<1 >+DE(2)+DE<3) 

DEAV=DEAV/3. 

DEIN1=DEAV*3. 

DDEXN2-<<DE< 1  >-DE<2>  > **2 . + < DE < 2 > -BE ( 3)  > **2 , + < DE < 3 > -DE < 
l>>**2.)/6. 

DEIN2=< <E< 1)-E<2) >**2 . + <E < 2 > -E < 3 > >**2 . + < E < 3 > -E < 1 > )**2. 

)  /4  » 

WRITE(3»175)SINl»SIN2fSIN3»DSIN2fDEINl»DDEIN2»DEIN2 
175  FORMAT ( 7G ) 

RETURN 

END 

SUBROUT INE ( FC  >  FY  r  ECT  »  PS1 >  PS2  >  PS3  > 

C  COM  C(60)»D<1»4> 

C<1 )=0.7 
C  <  2 ) =0  • 4/FC 

C( 3 >=1 400. *<F C/4430.  )**0,5 
C(4)=90.*3400./FC#*4. 

C ( 5 ) =0 . 043 
C<4)=0.4*3400/FC**2. 

C<7>=0. 15/FC**2. 

C ( 8 ) =0 . 05 

C<9)=15./FC**2,#(FC/3400)**1.5 

C<10>*1 .5E~3 

C< 1 1 )=1 « 2SE--4 

C<12)=0.2/FC 

C(13)=0.8/FC 

C  < 1 4 )  =  2  *  2E-5/F  C 

C ( 1 5 ) =25 ♦ 

C( 14)=1 .095 
C(17)=1.216 
C(18)=0.055 
C ( 19 ) =0  *  94 
C(20)=4300./FC**2. 

C ( 21 )=1 4 , 

C ( 22 ) = 1000 • 

C ( 23 ) =0 . 04 
C ( 24  >  =0 • 2#FC 
C(25)=9«1#F C/7020. 

C ( 24 ) =FC 
C<27)=2./FC 
C  <  28 )  =  3  •  E--3 
C  ( 2  9 )  =  0 «  5 
C  ( 3  0  )  =  2 
C(31 )®150. 

C  <  32 )  =  30 . 

C ( 33 ) =3500 . 

C  (  3  4 )  =  0 . 0  (3 
C< 35) =0.23 
C ( 34 ) = 1 . 5E-3 
C< 37 >=0.0125 
C(38)=3E-3 


r  r 


{  r 


177 

11523 

* 


C ( 39  >  =0  *  18 
C<40>=0.002 
C  <  4 1 )  =  1  • 05E-6 
C( 42) =0.001 
C  <  43 ) =  3E-  3 

C  (  44 )  =  4.  E6-I  (EC-4650)  *1000. 

C ( 45 ) =F Y 
C(46)= 2 9000000* 

C(47)=C<45>/C<46) 

C  <  48 ) =8 • *C ( 47 ) 

C<49)=80,*e<47> 

C  <  50  >=F'S1 
C  <  51  >  =F'S  2 
C  <  53 ) “PS3 

C<53>=5t*SQR(C<26) ) 

SOCRACK  SPACING 
AC=AREA  OF  SHEAR  PLANE 
CM-"MINIMUN  CONCRETE  COMER  TO  BAR 
NBL=NUMBER  OF  BAR  PER  LAYER 
NB=Nl)MBER  OF  LAYER  IN  SECTION 
BN=NET  BEAM  WIDTH 
C(54)=4 . 

DO  345  1=1*6 
D  < 1  *  I ) =0  ♦  0 
5  CONTINUE 
RETURN 
END 

SUBROUTINE  ROT  ATE  ( S .  S 1 » E  *  El  ,  HE  ) 

ZER0=0 . 0 
DO  176  1=1,4 
THE  < I ) =  ZERU 
>  CONTINUE 
ZER0= 1 

PREVIOUS  STRESSES  TO  PRINCIPAL  STRESSES 
CAL  PR  IN ( S 1  *  THE  *  ZERO ) 

CURRENT  STRESSES  TO  PRINCIPAL  STRESSES 
CALL  PRIN(S» THE, ZERO) 

C  PREVIOUS  STRAINS  TO  PRINCIPAL  STRAINS 
El <  4 ) =E 1 <  4 ) /2 . 

ZERO  =1, 

CALL  PR  IN (El, THE, ZERO) 

CURRENT  STRAIJNS  TO  PRINCIPAL  STRAINS 
E(4)=E(4)/2. 

CALL  PRIN(E,THE*ZERO) 

RETURN 

END 

SUBROUT  1 NE  FUNEND ( S , S 1 , E , E 1 , THE , DC , I STRE8) 
COM  C  <  60 ) , D ( 1 , 6 ) 

TOLET A  =  0  •  0 1 

TOLAM  =  TOLET  A 

ZERU=0  *  0 

DLAMP=ZERO 

DLAM=DLAMP 

DETA-ZERO 

DETAP  =  ZER(J 

I TE=  ZERO 

DO  177  1=1,4 

DS< I)*ZERO 

DE  <  I )  ~E  <  I )  --  f  t  <  I  ) 

23  ITE=ITE+1 


..  STRESSES 

STRESSES 

STRAINS 


f  t  ( I  ) 


-'OllJli.iliJ  nl  ,-j _ luJ  Aji/i.-t 


S  $DS< I ) 
e  I  <  i  > = e  1  ( t )  +he(  i  > 

CONTINUE 
THE  <  5 ) =  S 1(4) 
XLAHI=P< 1 > 1 ) IDE AM 
ET A I  =  D ! 1 » 3  >  4 DE T A 
XLAMPI=D<  1  f  2 )  IDLAHP 
ETAPI=D!  1  »4>+DElAP 


CALL  INVAR!  SI  »f.:  I » PE » S INI »  SIN2 » SI N3 » DSIN2»  DEIN1 » DDEIN2» 
PEIN2»ISTRES> 

ZER0=2  * 

DO  181  1  =  1.4 
181  STRESS! I >=BI ( I > 

CALL  PR I N<  S I RESS .  THE  f ZERO ) 

SMAX»8TRLSS!3> 

SHIN=S TRESS! 1 ) 

8MED"8TRES8(2> 

IF  !  ABS  (  SHAX  )  .IT.  1  E -4  )  SHAX=0 . 0 
IF  (  ABS !  SMED  )  ,t  T  .  IE -  4  )  SMEP=0 . 0 
IF! ABS ! SH I N ) .LT. I E -  4 ) SMtN=0 . 0 
G3“AI:*8  <  C  (  4 1  tSMIN )  *  .  93 - SttR  T  (  ME  I N2  ) 

IF  !G3.GT. 0,0)00  TO  11615 

XLPRI  =  SHIN*!ABS!G3)**!l/3)/U 4ABS!G3/C< 42) >**3, > ) 

GO  TO  11620 

11615  XLPRI=SNIN*<G3**!l/3)/<14AB3<63/C!42>  >**3.  )  ) 

11620  SLPR I =C ! 40  >  * ! 1 . -ABS ( XL AMP  I ) /C ( 43  )  ) 

Cl  =  l . /<  1  ,4C<31  )*XLAW1 ) 

E0=C ( 44 ) 

XK=C1*E0/!3.*!1 .  ~2.*C!39>  )  ) 

G=Cl#E0/!2.*! 1 ,4C!39)  )  ) 

XI  = < C < 50 ) *DEIN2/ ( C < 28 ) ¥*2 . 4DEIN2 > ) **2 . 
X1=X1+<XLAMI/C!38> >¥*2, 

Xl  =  Xl*(C!29)/( 1 .- C!27)*SIN1 )  ) 

Sl=l ,-XLAHI/C!33> 

H=C!25)*!SIN1/!C!26)-SIN1 ) >*»2. 

SH=1 .+ETAPI/C!34)E(ETAFI/C(35) )»#2. 

F5=C< 12>*SHIN¥! 1 .4C!13)*SMIN> 

F5=F5*!DEIN2*#0.25/ ! ABS!C< 14 ) *SH I N ) **0 . 254PE I N2**0 . 5 ) ) 


**3. 


<  24 )  ) 


F4=!DEIN2**0.25/!C!5)4PEIN2**0.25) >**3. 
G23=!DEIN2**0,25/!C!23)4DEIN2**0.5>  >**3. 
G22=14C!22)#!SHTN/!SHAX-C!24 ) > )**4 , 

G22=l/G22 

G21T0P=C!19)*!  SHED-SH IN )  /  (  SMED-C  (  24  )  )  - 1 . 

G21B0T=C!20)*! 1 .-C<21 > *ABS! SMIN ) / < SHAX-C ! 24 ) > )*!SHIN-C 


G21=!G21T0P/G21B0T)**1 .25 
PI =( SHED-SHIN)/! SHAX-C (24)  ) 

D2=C<  16>*tU#*!4./3.  )  - C(17) 

G1 1=C < 15 )*DEIN2*«0»: 5*01*02 

012=1 ,+!SHIN/(C! 1B)*!SHAX~C< 24) ) >  >**4. 

012=1 ,/G12 
G12=l , -012 
G1=G11*G12 
G2=G21*G22*G23 

F2T0P"C(3)*SQRT!DEIN2)*! 1 . 4ABS ! C ! 7 ) *SIN2 ) **0 . 254F5 > 
F2B0T»1  .-C!2)*SIN14ABS!C!9)*SIN2)**0.25*F4  C!4)*SIN3*D 
SIN2**. 125 

**!1.+G2> 

F2=F2T0P/F2B0T 
IF!SIN3,GT. 0.0)11850 

F1=C< 1)*! 1 .461 >/! 1 . 4C!6)*ABS!SIN3)**! 1 . /3 . ) #! 1 . 4G2 )  ) 

GO  TO  11860 

* 


/< 

r 

<> 

r 


rM 


•C<6>*A»S<SIN3>**<1  ./3.  >*(1.402)  > 
11870 


F3=l .4C<11)/<DEIN2*<1 .+C< 10)/ETAI**2. )  > 

SF  =  < 1 .4<C<32)*ETAI+C<33>*ETAI**2. ) /< 1 . +F2/C < 8 > ) >*F3 


11850  FI *C ( 1 ) * ( 1 . -G1 ) / < 1 
1 1860  F=F1 +F2 

IF (ETAI.GT.O.O)GO  TO 
F3  =  1 « 

GOTO  11880 

iils8  im 

C  CALCULATE  INCREMENTS  IN  LAMBDA, ETA, AND  INTRINSIC  TIME 
DPSI=SQRT ( DDEIN2 ) 

DPSIP  =  SQRT(pr:INl**2.) 

DETA~F*DPSt 
DETAP=H*DP8IP 
DCHI  =  DF.T  A/SF 
DCHXP«UETAP/SH 
DZ=DCHI/C<  36  > 
r.iZP  =  DCHIP/C<37) 

DLAMP=S1*X1*DPSI 

DELAMP=SLPRI*XLPRI*DPSI 

C  CALCULATE  DEVIATORIC  AND  VOLUMETRIC  STRESSES 
VSTRAI=DE(1 >+DE(2)+DE(3) 

VSTRES=SI<1 )+SI <2>4SI<3) 

IF (  iBS( VSTRES) .LT.  1 .E- 1 2 > VSTRES-0 . 0 
IF CABS ( VSTRA I ) > . LT . 1 . E-l 2 > VSTRA 1=0 . 0 
DO  191  1=1,3 

DS  IRAK  I  >=DE<  I) -VSTRA  1/3. 

DSTRES( I )=SI (I >-VSTRES/3. 

DSTRAI ( I )  =  DE ( 4  > 

DSTRES ( 4 ) =SI ( 4 ) 

DO  192  1=1,4 

B(I)=DSTRAI(I)/SQRT(4.*DDEIN2> 

XK0<I)=DSTRES(I)*F/(SF*C<36)  > 43 . XK* < S1*X1+SLPRI*XLPRI > 
IF(I.EQ.4.)XK0(I )= DSTRES ( I ) *F/ <  SF#C <  36 ) ) 
BNN=B<1)+B(2)+B<3) 

MAT  DC=ZER (4,4) 

DO  271  1=1,3 

DC(I, I)=XK+4,*G/3,-VSTRES*H/<3.*SH*C<37)  >-XKO(I)*B( I>- 
BNN/3 . ) 

DO  272  1=1,3 
K  =  1 4 1 

DO  273  J=K , 3 
DC< I , J)-XK-2. 

-BNN/3 . ) 

DC< J»I)=XK-2.*G/3.~VSTRES*H/<3, *SH*C<37) >-XNO< J)*<B(I> 
-BNN/3. > 

CONTINUE 
CONTINUE 
DO  274  J=1 , 3 

DC<  4*  J)=-XK0<4  )*(!<(  J) -BNN/3. > 

DC< J,4)~-XN0< J)*B<4) 

CONTINUE 

DC(4,4)=2.*6-XK0<4)*B<4) 

SOLVE  FOR  STRESSES  (1-PL  STRAIN?  2-PL  STRESS) 

IF ( ISTRES=0 ) GO  TO  12410 
GO  TO  (12360, 12410), ISTRES 
12360  DO  275  1=1,4 
DSTRES ( I ) =0 . 0 
DO  276  J  =  1 »  4 

DSTRES ( I )  =  DSTRES < I >  +DC ( I »  J  > *DE ( J ) 

276  CONTINUE 
275  CONTINUE 

GOTO  12500 
12410  DO  278  1=1,4 
DSTRES ( I) =0.0 


192 


271 


*G/3.-VSTRES*H/<3.*SH*C<37> )-XNO( I)*(B(J) 


2-~ 
272 


274 

C 


Zfjjb,  »****», 


i  * 


*DE<  J 

200 

12470 

2*2) 

C 

t: 

i 


12550 

I 


12614 


tf  <T  '?>nn  in  i:m;o 

LiU  200  J "  1 *4 

DSTRESS(I>=DSTREE(I >4!0C!I » J)-DC!I»2>*DC!2» J)/DC!2*2> ) 

CONTINUE 

CONTINUE 

HE  <  2  >  ( DC  ( 2  » 1 )  *DE  ( 1 )  +  DC  (  2  » 3  >  *DE  (  3 )  +  DC  ( 2  *  4 )  *DE  ( 4  > )  /DC  < 


DO  12550  1=1*6 
DS!I)=DSTRES5(I) 

CONTINUE 

IF ( ITE  *LT  « I )  60  TO  12614 
CDl.AM=ABS <  ( DL AM- PL AMI )  /DLAM) 

CBETA=A»S! < DETA-DETA1 > /DET A) 

IF ( CDLAH *  LT . TOLAM) »  AND • ( CDETA  *  LT . TOLETA )  DO  TO  12630 

IF! ITE  »0T  *7)  DO  TO  12630 

DLAM1=DLAM 

DETA1*DET  A 

60  TO  11523 


C 

C  CONVERT  STRESSES/STRAINS  TO  GLOBAL  COORDI 

NATES 

C 

IF ( NCR  » GT .0)G0T013500 
DO  12790  1=1*4 
E  <  I )  =E1  <  I  )+I»E<  I  > 

sim=s<  i) 

8(1 )*S( I >+DS( I ) 

1 2790  CONTINUE 
13500  RETURN 

SUBROUTINE  PRI N ( STRESS » THE  *  ZERO ) 

COMMON  C  <  60 ) 

IF ( ZERO  *  EG . 2 ♦ )  DO  TO  14120 
THE=STRESS< 1 ) -STRESS ( 3 ) 

IF (THE«NE*0 . )  60  TO  14060 
IF!THE*EQ*0.ANB.STRESS<4)  *0T  *0*  >  THE=PI/4. 

IF< THE. EG.O. AND. STRESS ( 4  >LT.0>  THE=3.#PI/4 . 

GO  TO  14080 

14060  THE=2.*STRESS!4)/THE 
TME*ATN<  THE )/2 . 

14080  STEMP=(STRESS(l)+STRESS(3))/2* 

TAUMAX=SQRT ( ( <  STRESS ( 1 )-STREE 

TAUMAX=SQRT(((STRESS!l)-STRESS!3>)/2*)*#2«+STRESS!4)#t 

2.  ) 

STRESS! 4  >  «0 • 

STRESS( 1 >=STEMP+TAUMAX 
STRESS! 3  >=S TEMP- TAUM AX 


IF ( ISTRES  *EG ♦ 0  >  60  TO  13500 
D!1*1)=D!1*1 HDLAM 
D! 1  *  2 ) =D ( 1  *  2 ) +DLAMP 
D!1*3)=D!1*3)+DETA 
P(1*4)=D(1,4  >  4DET  AP 


1 


14:  30 


14200 


14250 


1  4290 


14310 

14320 


14390 


14520 

14530 


14660 

C 

C 


,iU  ,u  1A,2,, 

8MAX=MIN1(STRESB(1)»STRESS<2) » STRESS  <  3  > ) 

IF  <  ZERO  » EG  *  1 ♦ )  RETURN 

H‘5n 

III“0 

no  14290  1=1*3 

IF  ( STRESS ( I ) , EG.SMAX)  GO  TO  14200 
IF  <STRESS<I) *  EG . SM IN  >  GO  TO  14250 
SME»=STRESS<  1  ) 

00  TO  14320 
11  =  11  +  1 

IF  <11.  EG. I)  IMED= I 
IF  (II.EQ.2)  GO  TO  14290 
IMAX= I 
GO  TO  14290 
111=111+1 

IF  (III, EG. 2)  IMED=I 

IF  (III. EG. 2)  GO  TO  14290 

TMJN=I 

CONTINUE 

IF ( 1 1 1 « EG . 3 ) .OR. (II .EG. 3) 

IF  (Ii>  ,NE,3,AND.II.NE.3>  GO  TO  14310 
1 MAX= I 
1 MED  =  I 
I M I N  =  I 

SMED=STRESS( IMED) 

STRESS ( 1 ) =SM I N 
STRESS(2)=SMEB 
STRESS! 3 ) =  SMAX 

Return 

SUBROUTINE  STEEL < ST , EPS i EPS1 * E r E 1 » THE r STRMAX > 

DIMENSION  DS ( 4  »  4  ) 

COMMON  C ( 60 ) »  D  < 1 »  6  > 

THIS  SUBROUTINE  IS  CALLED  BY  MATER1 . 

COMPUTE  STEEL  STIFFNESS  FOR  REPEATED  LOADING. 
DETERMINE  STEEL  MODULUS  IN  STEEL  LOCAL  DIRECTIONS  <R/ 

DO  14660  vl=  1  »  3 

I=J 

IF< (EPS(I >»EPSl(lf I)),LT,0.)  GO  TO  14520 
DES( J)=ABS(EPS( I ) > -ABS  <  EPS  1  ( 1  r  I  >  > 

GO  TO  14530 

DE8<J)*EPS(I)-EP3(1»I> 

IF  ( DES ( J ) « GT . 0 , )  GO  TO  14560 
ES(J)=C(46) 

GO  TO  14660 

STRAIN  HARDENING 
ES ( J ) =10000000 . 

STRMAX=EPS ( I ) 

CONTINUE 


DO  14700 

CALL  ZER ( DS  f  4 »  4 ) 

A  =  THE  <  2 ) 

DS<  1  »1)=C<50)*ES(1)*C0S(A)**4.+C(52>*ES(3)*SIN(A>**4 


DS(3»1>=(SIN<A)*C0S<A> )**2 . *( C (50 ) *ES ( 1 ) +C( 52 > *ES ( 3 ) > 


i 


1>^?,iM8IN(A)*C0S<A)  U(C(5°)*ES(  1  )*C0S(A)**2«-C<5 

1SIN(A)**2,> 

DS(2»2)=C<51 >*ES<2) 

DS<3»3>*C(50)*ES<1>*SIN<A>**4.+C<52)*ES<3>*C0S(A>**4. 
DS(4.4>»<2.*C0S(A>*SIN<A> > **2 . * < C <50>*ES < 1 ) +C ( 52 > *ES ( 3 

) ) 

DS<3» 4)=24*SIN(A>*C0S(A)*<C(50)*E8(l)*SIN(A)t*2.-C<52> 

*ES<3> 

1*C0S(A>**2.  > 

DS ( 3 . 1 ) *DS ( 1 » 3 ) 
nS(l»4)*DS(4»l> 

DS(4»3)*DS<3.4> 

Lit)  14830  1  =  1.4 
DE<X)"E<X)>E1<I> 

14830  CONTINUE 

DO  14880  1=1.4 
ST  < I ) *0 ♦ 

DO  14870  J  =  1 »  4 
ST(l>=ST<I>+DS<I» J)*DE(J> 

140/0  Continue 

14680  CONTINUE 
1 5500  RETURN 

SUBROUTINE  GLOBAL ( SX » A . B ) 

DIMENSION  DS(4.4) >80(4) 

CALL  Z E R ( D S  »  4 »  4  ) 

DS(1.1>=C0S<A>**2. 

DS<3»3)=C0S<A>**2. 

DS(t»3>=SIN(A>**2, 

PS( 3 » 1 ) «DS (1.3) 

DS(1»4)=-B*C0S(A)*SIN(A) 

D8<2»2)«1 . 

DS  ( 3  »  4 ) '-BS  < 1 »  4 ) 

DS ( 4 » 1 >  «C0S ( A ) #SIN  <  A ) 

IF(B.EQ. 1 .  >DS(4»1)=2.*C0S<A>*SIN<A> 

DS(4.3)=-DS<4. 1  ) 

DS(4.4)=C0S<A)*»2.-SIN(A)**2. 

DO  15660  1=1.4 
S0(I>*0. 

DO  15660  J=  1 » 4 
SO(I)*SO(I>+DS(I.J)*SX<J) 

15660  CONTINUE 

DO  15700  1=1.4 
SX(I)=SO(I) 

15700  CONTINUE 
RETURN 

SUBROUTINE  CRACHK ( S . S 1 . E . E 1 . CASTR » FT . NTCR » CRADIR ) 


15660 


15700 


CHECK  FOR  CRACKS 


THREE  DIRECTIONS 


DO  16080  1=1.3 

IF<CASTR(I).EQ.O.).AND.(S(I).LT.FT>  GO  TO  15830 
IF  (S(I).6T*FT> .AND. ( CASTR < I ) .EQ.O. )  GO  TO  15850 
IF  < CASTR ( I ) • NE . 0 . )  GO  TO  15900 
15830  CONTINUE 

C  NO  CRACKING  IN  THIS  DIRECTION 
GO  TO  16080 
15850  CONTINUE 

Prt-,  INITIAL  CRACKING  OCCURS  (TENSION  OR  COMPRFSION  STRAIN 
r  Ir.LP) 

PR0=(FT-S1<I))/<S<I)-S1(I)> 

CASTR< I ) =E1 < I ) +PRO#(S< I )~S1 ( I ) ) 

.  CRADIR*! 

* 


NTCR=NTCR+1 
CO  TO  16080 
15900  CONTINUE 

N  ELEMENT  PREVIOUSLY  CRACKED  IN  THIS  DIRECTION 

IF  <E1  <  I  >  .GT.CASTRT I) > .AND. <E< I) .GT.CASTR(I) >  00  TO  15 

950 

IF  <E1<I> .GT.CASTR(I) ) . AND . ( E < I ) . LT . CASTR ( I ) )  GO  TO  15 

980 

IF  <E1(I>.LT.CASTR<I>>.AND.<E<I>.LT.CASTR<I>>  GO  TO  16 

070 

IF  (Sl(I).LT.O. AND .S(I)«GT*0.)  GO  TO  16020 
15950  CONTINUE 

C  CRACK  REMAINS  OPEN 

NTCR=NTCR+1 
(JO  TO  16080 
15980  CONTINUE 
C  OPEN  CRACK  CLOSES 

NTCR=NTCR+1 
CASTR  < I )  =E  ( I ) 

GO  TO  16080 
16020  CONTINUE 
C  CLOSED  CRACK  OPENS 

P  R  0  =  - S 1 ( I )/(S< I >-Sl< I ) > 

CASTR ( 2 ) 

CASTR(I)=E1<I)+PR0*(E<I)-E1<I>) 

NTCR=NTCR+ 1 
GO  TO  16080 
16070  CONTINUE 

C  CLOSED  CRACK  REMAINS  CLOSED 

16080  CONTINUE 
16090  RETURN 

SUBROUTINE  ONECRA ( S » SI » E * El » SIG » SID1 * EPS » EPS 1 » 
*CASTR.C»TSTrTHE»DUN»DCRrA,FT, ICR> ISTRES) 

DIMENSION  ST<4)»ET(4)»DECR(4)»DEUN(J)»BE(4)»F(4»4)»SXX 
<  4 )  ,  DC  (  4 »  4  ) 

1 »  FD<  4 »  4 ) » THE  <  5 ) 

SC=C ( 54 ) 

AC=55 . 

I-ICR 

C 

T.TER  =  0 

IF  <  FT . EQ ♦ 0 ) GOTO 165 10 
PR0=(FT-S1(I))/(S(I)-S1(I)) 

C  UPDATE  STRESSES  TO  INCIPIENT  CRACKING 

DO  16420  I  =  1 »  4 
SXX ( I ) =EPS  < I > 

ST ( I ) KS1 ( I ) +PRO# ( S  < I > -SI ( I ) > 

SI ( I )  =ST  <  I ) 

El ( I )~E1 < I )  1PR04 < E ( I ) -El ( I )  > 

EPS  < I ) =EPS1 ( 1 . 1 ) +PRO#<  EPS< I ) -EPS1 < 1 1 1 ) ) 

16420  CONTINUE 
2  ♦ 

CALL ^GLOBAL ( ST  »A»B> 

DO  16500  1=1.4 
SIG ( I )  =ST ( I ) 

S ID1  <  1 . 1 )  =SI G ( I. ) 

IF ( I »NE» ICR)GO  TO  16480 
S(I)=S1<I) 

SI ( I )=0 

16480  EPS  1 ( 1 » I ) =EPS ( I > 

EPS  < I ) =SXX< I ) 

IF(I.EQ.2>EPS(2)=E1(2) 

16500  CONTINUE 


r 


r‘ 

r 


r 

.>*.• 

r 

r 


r 

r 


r 

r 

r 

r 


16560 

C 


16620 

16640 


16700 

C 

16800 


16860 

16880 


16920 

C 


16960 


DETERMINE.  PROF  ORT  ION  OP  101 AL  STRAIN  INCREMENT 
DO  16560  1=1*4 
DEUN< I)=(E<I>~El(I>)/2* 

DECR  I I ) *DEUN ( I ) 

DE<n*E<l)-El<X) 

IF! I *NE  *  2)G0  TO  16560 
DEON  1 2  >  =DE  <  2 ) 

DECR  <  2 ) *0 
CONTINUE 

DETERMINE  TOTAL  CRACKED  DIRECTIONS 
NCR  =  0 

DO  16620  1=1*3 

IF  <  CASTR (I).EG«0.0R.I.E0.2>  GO  TO  16620 
NCR=NCR  + 1. 

CONTINUE 
CONTINUE 
CALL  ZER ( F  *  4  *  4  > 

CALL  ZER ( DC  *  4  *  4 ) 

ITER  =ITER+1 

CALL  CRAST I I S 1  *  DCR  *  DECR  *  C . F  *  TST  *  THE » A »  NCR ) 

DUM=0 ♦ 

,  16700  1=1*4 
:T  ( I )  =DUN<  I )  +DEUN<  I  > 

ET1 < I )  =DUN ( I > 

CONTINUE 

CALL  FUNEND(S*S1.ET*ET1 *  THE , DC  * DUM *  NCR ) 

COMPUTE  UNCRACKED  STRAINS 
CALL  ZER < FD»  4  *  4 ) 

CALL  MATMU(F*DC»FD*4»  1*4) 

DO  16800  1=1*4 

FD( I  *  I )=FD< I  *  I ) +1 
CONTINUE 

CALL  INV I FD  *  DC »  4 ) 

DO  16880  1=1*4 

EoU16860°j=l *4 

DEUN( I)=DEUN<I)+DC(I* J)*DE( J> 

CONTINUE 

DECRI I )=DE< I )-DEUN( I ) 

CONTINUE 

IF  ( ITER *GT . 2 )  GO  TO  16920 
GO  TO  16650 
CONTINUE 

COMPUTE  STRESSES  IN  CRACKED  AND  UNCRACKED  CONCRETE 
ISTRES=DUM1 
DO  16960  1=1*4 
ET  < I ) =  DUN< I ) +DEUNI I ) 

ET1(I)=DUN(I) 

CONTINUE 

CALL  FUNEND(S*S1*ET»ET1 » THE » DC  * ISTRES » NCR > 

CALL  CRASTKS1  *  DCR  » DECR  *  C » F  *  TST  *  THE  *  A  *  NCR) 

CALL  ZER 1 FD * 4 ) 

CALL  INVIFD  *  DC  *  4 ) 

CALL  ZER  I SC  *  4 ) 

DO  1000  1=1*4 
DO  1000  J  = 1  *  4 
SCI  I  * J)=F< I  *  J)+FD< I  * J) 

NEXT  J 
NEXT  I 

CALL  ZER ( FD »  4  > 

CALL  INVIFD* SC *4) 


O 

r 


17030 


170A0 


8C0N( I > "0  » 

DO  17050  J=  1 » 4 

SCON  ( I  >  «SCON  (I  >  +  F  D  ( 1 » J )  *  DE.  <  J ) 

CONTINUE 

»cr<i>=dcr< i >+decr< i  > 

DUN<  I.  )«DUN<  I )  FDEUN(  I ) 

com  iNur: 

RETURN 

SUBROUTINE  CRAST I ( SI . E 1 . DECR , C » F * TST *  THE , A » NCR > 

1 1  F  -  0 

COMPUTE  FLEXIBILITY  COEFFICIENTS  FOR  CRACKED  PLANE 


18220 


C  COMPUTE  FLEXIBILITY  COEFFICIENTS  FOR  CRACKED  PLANE 

SL=C ( ) 

DB-C<56> 

AC»C(55) 

CM-C  ( 5 ) 

NB! =C ( 58 ) 

NB-C  <  59  > 

BN=C<40) 

FC=»C  <  26 ) 

FY=C(45) 

UO- 1.0 
FS=0  *  0 

CWsABS ( DCR  < 1 )  +  DECR ( 1 ) )*SC 
IF ( ITE ♦ EQ  *  1 )  CW  =  ABS(DCR(3)+DECR(3> >#SC 
18170  CD*ABS ( DECR ( 4 ) )*SC 
XKN==590000*C<50> 

XKN=590000*C<50> 

IF  <  I TE  *  CO  *  1 )  XKN=590000*Cv52> 

DUM~<3»4E5-XKN/CUMtl»  09E-7 
IF  (DUM.BT.O. )  60  TO  18220 
DUM=0. 

18220  CO-CW 

IF  (C0.LT.5E-3)  C0=,003 
AST=1000*AC/<3,9*(C0~ .002HDUM) 
g0Y^*92*DB**2.*SURr<FC+FY/lE6> 

0D0=DB*BN/NB1*< . 47+ , 54*CM/ ( BN/NB1**2 , +  DB ) > 

0DU=0DY 

IF  (VDO.LT.ODY)  gDU=0D0 
A1  ~2  ♦ 

A2=0  ♦ 

FS-TST ( 1 ) 

IF  (ITE, EG, 1>  FS=TST<3> 

FC=FC/1000 

Bl= 1 , -2 , ♦FS/FY 

IF  ( B1 «  LT  « 0  < )  B1 =0 , 

B2~F8#DB*SQRT ( FC ) / ( ,003*A1*FY> 
DST0P=-<A2*VDU/A1-2.*B2*CD))2.*B2 

DST0ps-(A2*gDU/Al-2.*B2*CD)*2 . *B2  +  2 , *B1*DB*0DU*SGRT ( FC 
)/( ,003*A1) 

IF  (CD.NE.O.)  60  TO  18350 
D8T*1000»*312.*NB1*DB**1.75*NB 
GO  TO  18380 

18350  DSB0T  =  (A2)|!VDU/A1-2,*B2*CD)**2,+4.*B1)KDB*SQRT(FC)*VDU*C 
D/< ,003*A1 ) 

DSB0T=2,*SGRT(BSB0T) 

DST=B2+DST0P/DSB0T 
DST=DST*1000,*NB 
F4-AST+US1 
F  4  - 1  #  /  F  4 


n  * 


I  * 


C*ABS  <  CW ) 

i f < r. . i  r .0.01  >t>o.ot 

F3*<3.9/XKN41.0VL-/*XNN/<i:W*(XKN*CW41,  )  )  >m/1000. 


Fglof* !i^6?CwS*<-. 63)4 < ,22*CU**< -.552)- 1 .034 >*FC 
F2B0T».1353*CW**<-.8>4< , 1 64*CWM<  -  .  707  >  - 1 . 379 > *FC 
F2«F2TQP*ALl/<F3BOT*XKN*tOOO. > 

F4»l ./<XKN*1000. > 

CALL  ZER<F .4.4) 

IF  (tTE.EQ.l)  GO  TO  20150 
F(1.1)«F1*AC/SC 
F(1.4)*F2*AC/SC 
F(4tl)*F3*AC/6C 
F<4»4)=F4*AC/SC 
IF  (NCR.EQ.l)  RETURN 
BB=C ( 5A ) 

AC»C<55> 

CH*2  *  625 
NB1*1 
NB  =  4 
BN*5 « 25 
FS~0 . 

ITE*ITE+1 
GO  TO  18170 


20150  F ( 3 » 3 ) “FI# AC /SC 
F (3.4  )  =  -F2*AC/8C 
F(4.3)«-F3*AC/SC 
F(4»4)eF(4»4) 4F4f AC/SC 
RETURN 

SUBROUTINE  ZER<WfN»M> 

DIMENSION  M(N»M> 

BO  1  1*1 fN 
BO  1  J®1 . M 
1  WU,J)=0.0 

RETURN 
END 

SUBROUTINE  MATMU(A»B»C»N»M»L) 
DIMENSION  A(N.M) »B(M.L)iC(N.L> 
DO  1  1=1 .N 
DO  1  J=1»L 
C(Ii J>*0.0 
DO  1  K  =  1 » M 

1  C<If J)=C(If J)4A<I»K)#B(Kf J) 


If-JURN 


